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A Gravitational Wave Background from Reheating after Hybrid Inflation 
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The reheating of the universe after hybrid inflation proceeds through the nucleation and subse- 
quent collision of large concentrations of energy density in the form of bubble-like structures moving 
at relativistic speeds. This generates a significant fraction of energy in the form of a stochastic 
background of gravitational waves, whose time evolution is determined by the successive stages of 
reheating: First, tachyonic preheating makes the amplitude of gravity waves grow exponentially 
fast. Second, bubble collisions add a new burst of gravitational radiation. Third, turbulent motions 
finally sets the end of gravitational waves production. From then on, these waves propagate unim- 
peded to us. We find that the fraction of energy density today in these primordial gravitational 
waves could be significant for GUT-scale models of inflation, although well beyond the frequency 
range sensitivity of gravitational wave observatories like LIGO, LISA or BBO. However, low-scale 
models could still produce a detectable signal at frequencies accessible to BBO or DECIGO. For 
comparison, we have also computed the analogous gravitational wave background from some chaotic 
inflation models and obtained results similar to those found by other groups. The discovery of such 
a background would open a new observational window into the very early universe, where the details 
of the process of reheating, i.e. the Big Bang, could be explored. Moreover, it could also serve in 
the future as a new experimental tool for testing the Inflationary Paradigm. 



I. INTRODUCTION 



Gravitational waves (GW) are ripples in space-time 
that travel at the speed of light, and whose emission by 
relativistic bodies represents a robust prediction of Gen- 
eral Relativity. The change in the orbital period of a 
binary pulsar known as PSR 1913-1-16 was used by Hulse 
and Taylor [l| to obtain indirect evidence of their ex- 
istence. Although gravitational radiation has not been 
directly detected yet, it is expected that the present uni- 
verse should be permeated by a diffuse background j)f 
GW of either an astrophysical or cosmological origin . 
Astrophysical sources, like the gravitational collapse of 
supernovae or the neutron star and black hole bina- 
ries' coalescence, produce a stochastic gravitational wave 
background (GWB) which can be understood as com- 
ing from unresolved point sources. On the other hand, 
among the backgrounds of cosmological origin, we find 
the approximately scale-invariant background produced 
during inflation Q , or the GWB generated at hypothet- 
ical early universe thermal phase transitions, from rel- 
ativistic motions of turbulent plasmas or from the de- 
cay of cosmic strings (^j ■ Fortunately, these backgrounds 
have very different spectral shapes and amplitudes that 
might, in the future, allow gravitational wave observa- 
tories like the Laser Interferometer Gravitational Wave 
Observatory (LIGO) the Laser Interferometer Space 
Antenna (LISA) the Big Bang Observer (BBO) Q 
or the Decihertz Interferometer Gravitational Wave Ob- 
servatory (DECIGO) 0], to disentangle their origin 0. 
Unfortunately, due to the weakness of gravity, this task 
will be extremely difficult, requiring a very high accuracy 
in order to distinguish one background from another. It is 
thus important to characterize as many different sources 
of GW as possible. 



There are, indeed, a series of constraints on some of 
these backgrounds, the most stringent one coming from 
the large-scale polarization anisotropics in the Cosmic 
Microwave Background (CMB), which may soon be mea- 
sured by Planck [8] , if the scale of inflation is sufficiently 
high. There are also constraints coming from Big Bang 
nucleosynthesis [9^], since such a background would con- 
tribute as a relativistic species to the expansion of the 
universe and thus increase the light element abundance. 
There is also a constraint coming from millisecond pul- 
sar timing [lo| . Furthermore, it has recently been pro- 
posed a new constraint on a GWB coming from CMB 
anisotropies [ll| . Most of these constraints come at very 
low frequencies, typically from 10~^* Hz to 10^* Hz, 
while present GW detectors work at frequencies of or- 
der 1-100 Hz, and planned observatories will ranee from 
10-3 Hz of LISA to 10^ Hz of Advanced-LIGO If 
early universe phenomena like flrst order phase transi- 
tions [12, [3 or cosmic turbulence [l3| occurred around 
the electro-weak (EW) scale, there is a chance than the 
GW detectors will measure the corresponding associated 
backgrounds. However, if such early universe processes 
occurred at the GUT scale, their corresponding back- 
grounds will go undetected by the actual detectors, since 
these cannot reach the required sensitivity in the high 
frequency range of 10^ — 10^ Hz, corresponding to the 
size of the causal horizon at that time. There are how- 
ever recent proposals to cover this range [T5| , which may 
become competitive in the near future. 

Moreover, present observations of the CMB 
anisotropies and the Large Scale Structure (LSS) 
distribution of matter seem to suggest that something 
like Inflation must have occurred in the very early uni- 
verse. We ignore what drove inflation and at what scale 
it took place. However, approximately scale-invariant 
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density perturbations, sourced by quantum fluctuations 
during inflation, seem to be the most satisfying expla- 
nation for the CMB anisotropies. Together with such 
scalar perturbations one also expects tensor perturba- 
tions (GW) to be produced, with an almost scale-free 
power spectra [3]. Because of the weakness of gravity, 
this primordial inflationary GWB should decouple from 
the rest of matter as soon as it is produced, and move 
freely through the Universe till today. At present, 
the biggest efforts employed in the search for these 
primordial GW come from the indirect effect that this 
background has on the B-mode polarization anisotropies 
of the CMB 0], rather than via direct detection. The 
detection of such a background is crucial for early 
universe cosmology because it would help to determine 
the absolute energy scale of inflation, a quantity that 
for the moment is still uncertain, and would open the 
exploration of physics at very high energies. 

In the early universe, after inflation, other backgrounds 
of GW could have been produced at shorter wavelengths, 
in a more 'classical' manner, rather than sourced by 
quantum fluctuations. In particular, whenever there are 
large and fast moving inhomogeneities in the matter dis- 
tribution, one expects the emission of GW. This is much 
like the situation in classical electrodynamics, but with 
some differences. At large distances from the source, the 
amplitude of the electromagnetic field Ai is expressed 
as the first derivative of the dipole moment di of the 
charge distribution of the source, Ai ~ di/cr, while the 
amplitude of the GW is given by the second deriva- 
tive of the quadrupolc moment of the mass distribution, 
hij ~ GQij/c^r. In both cases, the larger the velocity 
of the matter/charge distribution, the larger the ampli- 
tude of the radiation produced. Nevertheless, the main 
difference between the two cases is the weakness of the 
strength of gravity to that of electromagnetism. Thus, 
in order to produce a significant amount of gravitational 
radiation, it is required that the motion of huge masses 
occurs at speeds close to that of light for the case of astro- 
physical sources, or a very fast motion and high density 
contrasts in the continuous matter distribution for the 
case of cosmological sources. In fact, this is believed to 
be the situation at the end of inflaton, during the con- 
version of the huge energy density driving inflation into 
radiation and matter at the so-called reheating of the 
Universe Such an event corresponds to the actual 

Big Bang of the Standard Cosmological Model. 

Note that any background of GW coming from the 
early universe, if generated below Planck scale, immedi- 
ately decoupled upon production, as can be easily under- 
stood by the following dimensional analysis argument. 
Assuming that gravitons were in thermal equilibrium 
with the early universe plasma, at a temperature T, the 
gravitons' cross section should be of order a ~ G^T^ . 
Then, given the graviton number density n ^ and 
velocity v = 1, the gravitons' interaction rate should 
be r = (nav) - T^/M^. Since the Hubble rate is 
H - T'^/Mp, then V - H{T/Mpf, so gravitons could 



not be kept in equilibrium with the surrounding plasma 
for T < Mp. Therefore, GW produced well after Planck 
scale will always be decoupled from the plasma, and 
whatever their spectral signatures, they will retain their 
shape throughout the expansion of the Universe. Thus, 
the characteristic frequency and shape of the GWB gen- 
erated at a given time should contain information about 
the very early state of the Universe in which they were 
produced. Actually, it is conceivable that, in the not 
so far future, the detection of these GW backgrounds 
could be the only way we may have to infer the physi- 
cal conditions of the Universe at such high energy scales, 
which certainly no particle collider will ever reach. How- 
ever, the same reason that makes GW ideal probes of the 
early universe — the weakness of gravity — , is responsible 
for the extreme difficulties we have for their detection on 
Earth. For an extensive discussion see Ref. 

In a recent letter [3| we described the stochastic back- 
ground predicted to arise from reheating after hybrid in- 
flation. In this paper we study in detail the various pro- 
cesses involved in the production of such a background, 
whose detection could open a new window into the very 
early universe. In the future, this background could also 
serve as a new tool to discriminate among different in- 
flationary models, as each of these would give rise to a 
different GWB with very characteristic spectral features. 
The first stage of the energy conversion at the end of 
inflation, preheating [iGj . is known to be explosive and 
extremely violent, and quite often generates in less than 
a Hubble time the huge entropy measured today. The 
details of the dynamics of preheating depend very much 
on the model and are often very complicated because of 
the non-linear, non-perturbative and out-of-equilibrium 
character of the process itself. However, all the cases have 
in common that only specific resonance bands of the fields 
suffer an exponential instability, which makes their occu- 
pation numbers grow by many orders of magnitude. The 
shape and size of the spectral bands depend very much on 
the inflationary model. If one translates this picture into 
position-space, the highly populated modes correspond 
to large time-dependent inhomogeneities in the matter 
distributions which acts, in fact, as the source of GW we 
are looking for. 

For example, in single field chaotic inflation models, 
the coherent oscillations of the inflaton during preheat- 
ing generates, via parametric resonance, a population of 
highly occupied modes that behave like waves of mat- 
ter. They collide among themselves and their scattering 
leads to homogenization and local thermal equilibrium. 
These collisions occur in a highly relativistic and very 
asymmetric way, being responsible for the generation of 
a stochastic GWB [21, 22, 23] with a typical frequency 
today of the order of 10^ — 10^ Hz, corresponding to the 
present size of the causal horizon at the end of high-scale 
inflation. There is at present a couple of experiments 
searching for such a background, see Refs. fis'l, based of 
laser interferometry, as well as by resonant superconduct- 
ing microwave cavities [26| . 
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However, there are models like hybrid inflation in 
which the end of inflation is sudden [131 and the con- 
version into radiation occurs almost instantaneously. In- 
deed, since the work of Ref. [1^ we know that hybrid 
models preheat in an even more violent way than chaotic 
inflation models, via the spinodal instability of the sym- 
metry breaking field that triggers the end of inflation, 
irrespective of the couplings that this fleld may have to 
the rest of matter. Such a process is known as tachyonic 
preheating [H, [2^ and could be responsible for copious 
production of dark matter particles [30] , lepto and baryo- 
genesi s 13 Ij , topological defects ^] , primordial magnetic 
fields [34I, etc. 

It was speculated in Ref. [13] that in (low-scale) mod- 
els of hybrid inflation it might be possible to generate 
a stochastic GWB in the frequency range accessible to 
present detectors, if the scale of inflation was as low as 
^^inf ~ 1 TeV. However, the amplitude was estimated 
using the parametric resonance formalism of chaotic pre- 
heating, which may not be applicable in this case. In 
Ref. [29] (from now on referred to as paper I), it was 
shown that the process of symmetry breaking proceeds 
via the nucleation of dense bubble-like structures mov- 
ing at the speed of light, which collide and break up into 
smaller structures (see Figs. 7 and 8 of paper I). We con- 
jectured at that time that such collisions would be a very 
strong source of GW, analogous to the GW production 
associated with strongly first order phase transitions [13] ■ 
As we will show in this paper, this is indeed the case dur- 
ing the nucleation, collision and subsequent rescattering 
of the initial bubble- like structures produced after hybrid 
inflation. During the different stages of reheating in this 
model, gravity waves are generated and ampliflcd until 
the Universe finally thermalizes and enters into the initial 
radiation era of the Standard Model of Cosmology. From 
that moment until now, during the whole thermal history 
of the expansion of the universe, this cosmic GWB will 
be redshifted as a radiation-like fluid, totally decoupled 
from any other energy-matter content of the universe, 
such that today's ratio of energy stored in these GW 
to that in radiation, could range from U^^h^ ~ 10"^, 
peacked around / ~ 10^ Hz for the high-scale models, 
to il^^h'^ ~ 10~^^, peacked around / ~ 1 Hz for the 
low-scale models. 

Finally, let us mention that since the first paper by 
Khlebnikov and Tkachev |2l|, studing the GWB pro- 
duced at reheating after chaotic infiation, it seems ap- 
propriate to reanalyze this topic in a more detailed way. 
The idea was extended to hybrid inflation in [2^ . but 
within the parametric resonance formalism. It was also 
revisited very recently in Ref. [H, ^] for the and 
m^(/)^ chaotic scenarios, and reanalysed again for hybrid 
inflation in Ref. [isj . this time using the new formalism 
of tachyonic preheating [H, [13] . Because of the increase 
in computer power of the last few years, we are now able 
to perform precise simulations of the reheating process 
in a reasonable time scale. Moreover, understanding of 
reheating has improved, while gravitational waves detec- 



tors are beginning to attain the aimed sensitivity 4] . Fur- 
thermore, since these cosmic GWBs could serve as a deep 
probe into the very early universe, we should character- 
ize in the most detailed way the information that we will 
be able to extract from them. 

The paper is divided as follows. In Section II we briefly 
review the hybrid model of infiation. Section HI is dedi- 
cated to our approach for extracting the power spectrum 
of GW from reheating. In section IV, we give a detailed 
account of the lattice simulations performed with two 
codes: our own FORTRAN parallelized computer code 
(running in the MareNostrum supercomputer (33j and 
in our UAM-IFT cluster [13]), as well as with a modi- 
fied version of the publicly available C-I--I- package LAT- 
TICEEASY m . Section V is dedicated to study the spa- 
tial distribution of the production of gravitational waves. 
In Section VI, we reproduce as a crosscheck, some of the 
results of [U, [13, [IB] concerning the GWB produced at 
reheating after chaotic inflation models. Finally, in sec- 
tion VII, we give our conclusions and perspectives for the 
future. 

II. THE HYBRID MODEL 

Hybrid infiation models [13] arise in theories of par- 
ticle physics with symmetry breaking fields ('Higgses') 
coupled to fiat directions, and are present in many ex- 
tensions of the Standard Model, _both in string theory 
and in supersymmetric theories [l^. The potential in 
these models is given by 

t/(a>,x) = A(^<i>t<i>-^) +9\'<p^<f + l^i\\ (1) 

where the contraction (f>1^$ should be understood as the 
trace Tr$^$ = ^](^]^. Inflation occurs along the lifted 
flat direction, satisfying the slow-roll conditions thanks to 
a large vacuum energy po = A?;^/4. Inflation ends when 
the inflaton x falls below a critical value and the sym- 
metry breaking field acquires a negative mass squared, 
which triggers the breaking of the symmetry and ends in 
the true vacuum, ]</)] — v, within a Hubble time. These 
models do not require small couplings in order to gener- 
ate the observed CMB anisotropics; e.g. a working model 
with GUT scale symmetry breaking, v — 10~^ Mp, with 
a Higgs self-coupling A and a Higgs-infiaton coupling 
g given by g = •\/2A = 0.05, satisfies all CMB con- 
straints [2u\, and predicts a tiny tensor contribution to 
the CMB polarization. The main advantage of hybrid 
models is that, while most chaotic infiation models can 
only occur at high scales, with Planck scale values for the 
inflaton, and ^j^^'' - lO^^ GeV, one can choose the scale 
of infiation in hybrid models to range from GUT scales 
all the way down to GeV scales, while the Higgs v.e.v. 
can range from Planck scale, v = Mp, to the Electroweak 
scale, V = 246 GeV, see Ref. [13, Hj- 

There are a series of constraints that a hybrid infla- 
tion model should satisfy in order to be in agreement 
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with observations. First of all, inflation should end in 
less than one e-fold, otherwise unacceptable black holes 



would form 



This can be written as the waterfall 



condition [27|, XmMp :$> M"^, which becomes 



m V 



(2) 



Then there is the condition, known as the COBE normal- 
ization, that the scalar amplitude should satisfy As = 
/2TTij) ~ 5 X 10~^, which gets translated into 



5 = (n - 1) 



Mf 



X io-4e("-i)^/2 



(3) 



(4) 



as well as the spectral tilt, 

1 Ml 
n-l = - — ^<0.05 

and finally the fact that we have not seen so far any 
tensor (gravitational wave) contribution in the CMB 
anisotropics, r — A^/A^ < 0.3, imposes the constraint 



AV4<2x10-=^^ 



(5) 



Taking all these conditions together, we find that a 
model with v = O.lA/p is probably ruled out, while one 
with V — 0.01 Mp is perfectly consistent with all ob- 
servations, and with reasonable values of the coupling 
constants, e.g. g — 4 x 10^'' and A = 10^'^. How- 
ever, the lower is the scale of infiation, the more diffi- 
cult it is to accommodate the amplitude of the CMB 
anisotropics with reasonable values of the parameters. 
For a scale of inflation as low as 10^^ GeV, one must 
signiflcantly finetune the couplings, although there are 
low scale models based on supersymmetric extensions of 
the standard model which can provide a good match to 
observations [svj . 

In the following sections we will show how efficient 
is the production of GW at reheating after hybrid in- 
flation, using both analytical estimates and numerical 
simulations to derive the amplitude of the present day 
GWB. Reheating in hybrid inflation goes through 
four well defined regimes: first, the exponential growth 
of long wave modes of the Higgs field via spinodal in- 
stability, which drives the explosive growth of all parti- 
cles coupled to it, from scalars [SsJ to gauge fields [3l| 
and fermions [30]; second, the nucleation and collision 
of high density contrast and highly relativistic bubble- 
like structures associated with the peaks of a Gaussian 
random field like the Higgs, see paper I; third, the tur- 
bulent regime that ensues after all these 'bubbles' have 
collided and the energy density in all fields cascades to- 
wards high momentum modes; finally, thermalization of 
all modes when local thermal and chemical equilibrium 
induces equipartition. The first three stages can be stud- 
ied in detailed lattice simulations thanks to the semi- 
classical character of the process of preheating [ll] , while 
the last stage is intrinsically quantum and has never been 
studied in the lattice. 



III. GRAVITATIONAL WAVE PRODUCTION 

Our main purpose in this paper is to study the details 
of the stochastic GWB produced during the reheating of 
the universe after hybrid inflation (sections IH, IV and 
V). However, we also study, albeit very briefly, the anal- 
ogous background from reheating in some simple chaotic 
models (section VI). Thus, in this section we derive a 
general formalism for extracting the GW power spectrum 
in any scenario of reheating within the (flat) Friedman- 
Robertson- Walker (FRW) universe. The formalism will 
be simplifled when applied to scenarios in which we can 
neglect the expansion of the universe, like in the case of 
most Hybrid models. 

A theory with an inflaton scalar field x interacting with 
other Bose fields (pa, can be described by 



R 



WnG 



Vi(t>,x) (6) 



with R the Ricci scalar. For hybrid models, we consider 
a generic symmetry breaking 'Higgs' field with Nc real 
components. Thus, we can take — ^ J2a ^1 = 
in ([T]), with a running for the number of Higgs' compo- 
nents, e.g. Nc = I for a real scalar Higgs, Nc = 2 for a 
complex scalar Higgs or Nc — 4 for a SU{2) Higgs, etc. 
The effective potential ^ then becomes 



v'y + g\M' + lt^'x'. (7) 



For chaotic scenarios, we consider a massless scalar field 
interacting with the inflaton via 



(8) 



with V{x) the inflaton's potential. Concerning the sim- 
ulations we show in this paper, we concentrate in the 
Nc — 4: case for the hybrid model and consider a poten- 
tial V{x) = jX^ for the chaotic scenario. 

The classical equations of motion of the inflaton and 
the other Bose fields are 



X + 3i?X - \ 



1 



dV_ 

dx 







dV 







(9) 
(10) 



with H = d/a. 

Gravitational Waves are represented here by a 
transverse-traceless (TT) gauge-invariant metric pertur- 
bation, hij , on top of the flat FRW space 



-dt^ + a^{t) {Sij + h^j) dx^dx^ 



(11) 



with a{t) the scale factor and the tensor perturbations 
verifying dihij — ha = 0. In the following, we will raise 
or low indices of the metric perturbations with the delta 



Kronecker 6ij, so hij 



h^^ and so on. The Ein- 



stein field equations can be splitted into the background 
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(0) 



8ttGT^°J and the perturbed SG^u = 8,ttG6T. 



G 

equations. The background equations describe the evo- 
lution of the flat FRW universe through 



with diFi 



H 
4^ 



(13) 

where any term in the r.h.s. of (fT2|) and (I13|l . should be 
understood as spatially averaged. 

On the other hand, the perturbed Einstein equations 
describe the evolution of the tensor perturbations ^SQj as 



IQirGWi 



(14) 



with 9^11^^ = Ilii = 0. The source of the GW, Ilij, con- 
tributed by both the inflaton and the other scalar fields, 
will be just the transverse-traceless part of the (spatial- 
spatial) components of the total anisotropic stress-tensor 

^^^J.u = ^ [^^J,X^uX + d^,(j)ady(t)a + g,,,yi£' - (p))] , (15) 

where £(%, 0a) is the lagrangian ^ and (p) is the back- 
ground homogeneous pressure. As we will explain in the 
next subsection, when extracting the TT part of (|15p . 
the term proportional to (7^1/ in the r.h.s of (|15p. will be 
dropped out from the GW equations of motion. Thus, 
the effective source of the GW will be just given by the 
TT part of the gradient terms d^x^^vX + dti4>adu4>a- 

In summary, Eqs. ([5|)- PII)) . together with Eqs. IT^ - 
(I13p . describe the coupled dynamics of reheating in any 
inflationary scenario, while Eq. ((T4)) describe the pro- 
duction of GW in each of those scenarios. In this pa- 
per we use lattice simulations to study the generation 
of GW during reheating after inflation. Specific details 
on this are given in section IV, but let us just mention 
here that our approach is to solve the evolution of the 
gravitational waves simultaneously with the dynamics of 
the scalar fields, in a discretized lattice with periodic 
boundary conditions. We assume initial quantum fluc- 
tuations for all flelds and only a zero mode for the infla- 
ton. Moreover, we also included the GW backreaction 
on the scalar fields' evolution via the gradient terms, 
h^^y iX^ jX + h^''^i4'a^j4'a and confirmed that, for all 
practical purposes, these are negligible throughout GW 
production. 



A. The Transverse-Traceless Gauge 



A generic (spatial-spatial) metric perturbation Shij has 
six independent degrees of freedom, whose contributions 
can be split into scalar, vector and tensor metric pertur- 
bations 13911 



and dihij 



hi, 



0. By choosing a 



transverse-traceless stress-tensor source 11^ j , we can elim- 
inate all the degrees of freedom (d.o.f.) but the pure TT 
part, hij, which represent the only physical d.o.f which 
propagate and carry energy out of the source (GW). If 
we had chosen only a traceless but non-transverse stress 
source, we could have eliminated the scalar d.o.f. ip and 
absorbed E into the scalar field perturbation, but we 
would still be left with a vector field Fi also sourced by 
the (traceless but non-transverse) anisotropic stress ten- 
sor, thus giving rise to a vorticity divergence-less field 
Vi. However, since the initial conditions are those of a 
scalar Gaussian random field (see section IV), even in 
that case of a non-transverse but traceless stress source, 
the mean vorticity of the subsequent matter distribution, 
averaged over a sufliciently large volume, should be zero 
(although locally we do have vortices of the Higgs field, 
see Refs. [3l|,[33|), since vortices with opposite chirality 
cancel eachother. This means that in such a case, al- 
though d^Ilij ^ 0, and thus d^Shij ^ 0, we could still re- 
cover the TT component when averaging over the whole 
realization. 

For practical purposes, we will consider from the beg- 
gining the TT part of the anisotropic stress-tensor, ensur- 
ing this way that we only source the physical d.o.f. that 
represent GW. The equations of motion of the TT metric 
perturbations are then given by Eq. Note that for 

a non-transverse source the equations would have been 
much more complicated, so the advantage of using the 
TT part from the beginning is clear. The disadvantage 
arises because obtaining the TT part of a tensor in conflg- 
uration space is very demanding in computational time. 
However, as we explain next, we will use a method by 
which we can circunvent this issue. 

Let us switch to Fourier space. Using the convention 



/(k) 



1 



(27r)3/2 



+zkx 



/(x) 



(17) 



the GW equations (fTi|) in Fourier space read 



/ly {t, k) + 3Hh,j (i, k) + — /ly it, k) = IGTTG Hy (t , k) , 

(18) 

where k — |k|. Assuming no GW at the beginnig of re- 
heating (i.e. the end of inflation te), the initial conditions 
are hij{te) — hij{te) — 0, so the solution to Eq. psp for 
t > tf, will be just given by a causal convolution with an 
appropriate green function G{t,t'), 



hij{t,k) = ISttG^ rft'G(t,t')n,y(i',k) 



(19) 



F, 



(16) 



Therefore, all we need to know for computing the GW is 
the TT part of the stress-tensor, Hy, and the Green func- 
tion G{t' , t). However, as we will demonstrate shortly, we 
have used a numerical method by which we don't even 
need to know the actual form of G{t' , t). To see this, let 
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us extract the TT part of the total stress-tensor. Given 
the symmetric anisotropic stress-tensor T^i^ (jl5|) , we can 
easily obtain the TT part of its spatial components in 
momentum space, ny(k). Using the spatial projection 
operators Pij = 5ij — kikj, with ki = ki/k, then (40| 



where 



A,,,i,„(k) = (F,,(k)P,„(k) - iPi,(k)Pz„(k) 



(20) 



(21) 



Thus, one can easily see that, at any time t, kiHij (k, t) = 
n*(k, t) — 0, as required, thanks to the identities Pijkj = 
and P^jPjjYi — Pi7n- 

Note that the TT tensor, Ily, is just a linear combi- 
nation of the components of non-traceless nor-transverse 
tensor T^ (fTS)) . while the solution (fT9)) is just hnear in 
Hij . Therefore, we can write the TT tensor perturbations 
(i.e. the gravitational waves) as 

h,j (i , k) = Ay- (k)u„- (i , k) , (22) 

with Uij{t, k) the Fourier transform of the solution of the 
following equation 



1 2 

ili-j + SHuij Ui 



167rGT, 



(23) 



This Eq. is nothing but Eq. HH), sourced with the 
complete Ty (fT5|) . instead of with its TT part, Ily. Of 
course, Eq. ((23|) contains unphysical (gauge) d.o.f.; how- 
ever, in order to obtain the real physical TT d.o.f., hij , we 
can evolve Eq. (|23p in configuration space, Fourier trans- 
form its solution and apply the projector ([2T|) as in ([22l) . 
This way we can obtain in momentum space, at any mo- 
ment of the evolution, the physical TT d.o.f. that repre- 
sent GW, hij. Whenever needed, we can Fourier trans- 
form back to configuration space and obtain the spatial 
distribution of the gravitational waves. 

Moreover, since the second term of the r.h.s of the 
total stress-tensor Ty is proportional to gij = Sij + hij , 
see ([TS]) . when aplying the TT projector ((2T|) . the part 
with the Sij just drops out, simply because it is a pure 
trace, while the other part contributes with a term — {C — 
{p))h.ij in the l.h.s of Eq.([TH]). However, (£ — (p)) is of 
the same order as the metric perturbation ~ 0{h), so this 
extra term is second order in the gravitational coupling 
and it can be neglected in the GW Eqs. p^ . This way, 
the effective source in Eq. ([25)1 is just the gradient terms 
of both the inflaton and the other scalar fields, 



T,; 



(V«xVjX + V,(/)aVj</.a) 



(24) 



Therefore, the effective source of the physical GW, will 
be just the TT part of as we had already mentioned 
before. 

Alternatively, one could evolve the equation of the TT 
tensor perturbation in configuration space, Eq. (jl4[) . with 
the source given by 



nij(x,t) 



1 



(27r)3/2 



d3ke-*''A,,,,,„(k)Ti„(k,i), (25) 



such that diliij{'x.,t) — nii(x, t) at any time. So, at 
each time step of the evolution of the fields, one would 
have first to compute (the gradient part of) T;„i ([M]) in 
configuration space, then Fourier transform it to momen- 
tum space, substitute in Eq. (^5]) and perform the inte- 
gral. However, proceeding as we suggested above, there 
is no need to perform the integral, nor Fourier transform 
the fields at each time step, but rather only at those times 
at which we want to measure the GW spectrum. The vi- 
ability of our method relies in the following observation. 
To compute the GW we could, first of all, project the 
TT part of the source (PS)) . and second, solve Eq. pi)) . 
However, we achieve the same result if we commute these 
two operations such that, first, we solve Eq. and sec- 
ond, we apply the TT projector to the solution (P^ . We 
have found this commuting procedure very useful, since 
we are able to extract the spectra or the spatial distri- 
bution of the GW at any desired time, saving a great 
amount of computing time. Most importantly, with this 
procedure we can take into account backreaction simul- 
taneously with the fields evolution. 

In summary, for solving the dynamics of reheating of 
a particular inflationary model, we evolve Eqs. ([5|)- (fTU)) 
in the lattice, together with Eqs. (fT^ - P^ . while for 
the GWs we solve Eq. (|23| . Then, only when required, 
we Fourier transform the solution of Eq. (pS]) and then 
apply in order to recover the physical transverse- 
traceless d.o.f representing the GW. From there, one can 
easily build the GW spectra or take a snapshot of spatial 
distribution of the gravitational waves. 



B. The energy density in GW 

The energy-momentum tensor of the GW is given 
by 



1 



327rG 



(26) 



where hij are the TT tensor perturbations solution of 
Eg. p^ . The expectation value (...)v is taken over a 
region of sufficiently large volume V = to encompass 
enough physical curvature to have a gauge-invariant mea- 
sure of the GW energy-momentum tensor. 

The GW energy density will be just p^^ — <oo, so 



Pgw 
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d^x hij {t, x.)hij {t, x) 
d3k/i,j(t,k)/i^(t,k), (27) 



where in the last step we Fourier transformed each 
hij and used the integral definition of the Dirac delta 
(27r)35(3)(k) =Jd^K e-'kx^ 

We can always write the scalar product in (|27p in terms 
of the (Fourier transformed) solution u;„i of the Ea. (P5)) . 
by just using the spatial projectors ((2T|) 



hah. 



A 
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(28) 
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where we have used the fact that f^ij^iml^im.rs — ^ij.rs- 
This way, we can express the GW energy density as 



1 



PGW 
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k^dk I dr!Ay-i™(k)uy(t,k)w*„(i,k). (29) 



From here, we can also compute the power spectrum per 
logarithmic frequency interval in GW, normalized to the 
critical density pc, as 



where 



^gw (^) 



]_dpgw_ 

Pc d logfc 
32nGL3pc 



(/) 



(30) 



df! Ay- ,„(k)M,,(i, k)ur,„(t, k) (31) 



We have checked explicitely in the simulations that the 
argument of the angular integral of ([31]) is independent 
of the directions in k-space. Thus, whenever we plot 
the GW spectrum of any model, we will be showing the 
amplitude of the spectrum (per each mode k) as obtained 
after avaraging over all the the directions in momentum 
space, 



^gw (^) 



fc3 



SGL^pc 



A,,,i™(k)u,,(t,k)uL(t,k)) (32) 

47r 



where {f),^^^Jfdn. 

Finally, we must address the fact that the frequency 
range, for a GWB produced in the early universe, will be 
redshifted today. We should calculate the characteristic 
physical wavenumber of the present GW spectrum, which 
is redshifted from any time t during GW production. 
This is a key point, since a relatively long period of turbu- 
lence will develop after preheating, which could change 
the amplitude of the GWB and shift the frequency at 
which the spectra peaks. So let us distinguish four char- 
acteristic times: the end of inflation, tg] the time when 
GW production stops; the time tr when the universe fi- 
nally reheats and enters into the radiation era; and today, 
to.^ Thus, today's frequency /o is related to the physi- 
cal wavenumber kt at any time t of GW production, via 
fo = (at/ao)(fct/27r), with oq and at, the scale factor to- 
day and at the time respectively. Thermal equilibrium 
was established at some temperature Tr, at time tr > t. 
The Hubble rate at that time was MpHr = (87r/3)pr, 
with pr = grTr'^Tr /SO the relativistic energy density and 
gr the effective number of relativistic degrees of freedom 



^ Note, however, that after thermahzation there is still a small 
production of GW from the thermal plasma, but this can be 
ignored for all practical purposes. 



at temperature Tr- Since then, the scale factor has in- 
creased as flr/ao = {go^s/gr,sy^^{TQ/Tr), with gi^s the 
effective entropic degrees of freedom at time U, and To 
today's CMB temperature. Putting all together, 



fSlT^gr 

\ 90 



go,s 

gr,s 



Tn 



-) f , (33) 
ar I Zn 



where we have used the fact that the physical wave num- 
ber kt at any time t during GW production, is related to 
the comoving wavenumber k through kt — {ae/at)k with 
the normalization = I. 

From now on, we will be concerned with hybrid infla- 
tion, leaving chaotic inflation for section VI. Within the 
hybrid scenario, we will analyse the dependence of the 
shape and amplitude of the produced GWB on the scale 
of hybrid inflation, and more specifically on the v.e.v. 
of the Higgs field triggering the end of inflation. The 
initial conditions are carefully treated following the pre- 
scription adopted in paper I, as explained in section IV. 
Given the natural frequency at hand in hybrid models, 
m — \/Xv, whose inverse sets the characteristic time 
scale during the first stages of reheating, it happens that 
as long as w <C Mp, the Hubble rate H ^ ^/X{v'^ / Mp) 
is much smaller than such a frequency, H <^ m. Indeed, 
all the initial vacuum energy pq gets typically converted 
into radiation in less than a Hubble time, in just a few 
time steps. Therefore, we should be able to ignore 
the dilution due to the expansion of the universe during 
the production of GW, at least during the first stages 
of reheating. However, as we will see, the turbulent be- 
haviour developed after those first stages, could last for 
much longer than an e-fold, in which case we will have 
to take into account the expansion of the universe. Our 
approach will be first to ignore the expansion of the Uni- 
verse and later see how we can account for corrections 
if needed. Thus, we set the scale factor a — 1 and the 
Hubble rate H — and H — 0. As we will see later, our 
approach of neglecting the expansion for the time of GW 
production, will be completely justified a posteriori. 

The coupled evolution equations that we have to solve 
numerically in a lattice for the hybrid model are 



IGTrGT; 



x-^\+{g^\^f + f^^)x^o 



(34) 
(35) 
(36) 



with Ty given by Eq. ([24ll with the scale factor a = 1 . We 
have explicitly checked in our computer simulations that 
the backreaction of the gravity waves into the dynamics 
of both the inflaton and the Higgs fields is negligible and 
can be safely ignored. We thus omit the backreaction 
terms in the above equations. 

We evaluate during the evolution of the system the 
mean field values, as well as the different energy compo- 
nents. As shown in Fig.[Tl the Higgs field grows towards 
the true vacuum and the inflaton moves towards the min- 
imum of its potential and oscillates around it. We have 
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FIG. 1: Time evolution of the mean field values of the Higgs 
and the Inflaton, the former normalized to its v. e.v., the latter 
normalized to its critical value xo = m/g- This is just a 
specific realization with A*' — 128, Pmin = 0.1m, a = 0.48m~^, 
V = IQ-^Mp and 5^ = 2A = 0.25. 

checked that the sum of the averaged gradient, kinetic 
and potential energies (contributed by both the inflaton 
and the Higgs), remains constant during reheating, as ex- 
pected, since the expansion of the universe is irrelevant 
in this model. We have also checked that the time evo- 
lution of the different energy components is the same for 
different lattices, i.e. changing the number of points N 
of the lattice, of the minimum momentum Pmin — in/L 
or of the lattice spacing a ~ L/N, with L the lattice 
size, as long as the product ma < 0.5; for a detailed 
discussion of lattice scales see paper I. Looking at the 
time evolution of the Higgs' v. e.v. in Fig. [U three stages 
can be distinguished. First, an exponential growth of the 
v.e.v. towards the true vacuum. This is driven by the 
tachyonic instability of the long- wave modes of the Higgs 
field, that makes the spatial distribution of this field to 
form lumps and bubble-fike structures [2l,[2§|. Second, 
the Higgs field oscillates around the true vacuum, as the 
Higgs' bubbles collide and scatter off eachother. Third, 
a period of turbulence is reached, during which the in- 
fiaton oscillates around its minimum and the Higgs sits 
in the true vacuum. For a detailed description of the 
dynamics of these fields see Ref. [1^. Here we will be 
only concerned with the details of the gravitational wave 
production. 

The initial energy density at the end of hybrid inflation 
is given by po = m^?;^/4, with m? = Aw^, so the fractional 
energy density in gravitational waves is 



4t 



00 



pa rn? 



where ( hij h}^ 



o n2 2 (^'J^" 



(37) 



defined as a volume average like 



Y / (Pxhij h^^ , is extracted from the simulations as 



--y I dlogfc/c3^Ay to,(k)u„ (t,k)Mr„(t,k)^^^ (38) 
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FIG. 2: The time evolution of the different types of energy (ki- 
netic, gradient, potential, anisotropic components and grav- 
itational waves for different lattices), normalized to the ini- 
tial vacuum energy, after hybrid inffation, for a model with 
V = 10""^ Mp. One can clearly distinguish here three stages: 
tachyonic growth, bubble collisions and turbulence. 



where My (i, k) is the Fourier transform of the solution 
of Eq. (|34p . Then, we can compute the corresponding 
density parameter today (with fij-ad — 3.5 x 10^^) 



O /i2 



2Gv^rn?V 



dlogk ( Ayy,„(k)uy (t, k)ul„^{t, k) ) 



/ 47r 



(39) 



which has assumed that all the vacuum energy po gets 
converted into radiation, an approximation which is al- 
ways valid in generic hybrid inflation models with v <ti 
Mp, and thus H m — ^/Xv. 

We have shown in Fig. [2] the evolution in time of the 
fraction of energy density in GW. The first (tachyonic) 
stage is clearly visible, with a (logarithmic) slope twice 
that of the anisotropic tensor Hy- . Then there is a small 
plateau corresponding to the production of GW from 
bubble collisions; and finally there is the slow growth 
due to turbulence. In the next section we will describe 
in detail the most significant features appearing at each 
stage. 

Note that in the case that H ^ m, the maximal pro- 
duction of GW occurs in less than a Hubble time, soon 
after symmetry breaking, while turbulence lasts several 
decades in time units of m^^. Therefore, we can safely 
ignore the dilution due to the Hubble expansion, up to 
times much greater than those of the tachyonic instabil- 
ity. Eventually the universe reheats and the energy in 
gravitational waves redshifts like radiation thereafter. 

To compute the power spectrum per logarithmic fre- 
quency interval in GW, rjgu,(/), we just have to use ([3T|). 
Moreover, since gravitational waves below Planck scale 
remain decoupled from the plasma immediately after pro- 
duction, we can evaluate the power spectrum today from 
that obtained at reheating by converting the wavenum- 
ber k into frequency /. Simply using Eq. p3p . with 
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FIG. 3: We show here the comparison between the power 
spectrum of gravitational waves obtained with increasing lat- 
tice resolution, to prove the robustness of our method. The 
different realizations are characterized by the the minimum 
lattice momentum (pmin) and the lattice spacing (ma). The 
growth is shown in steps of mAt = 1 up to mt — 30, and then 
in and mAt = 5 steps up to mt — 60. 

9r,s/9o,s ~ 100, gr.s ~ 5r and Og ~ a*, then 

/ - 6 X IQi" Hz - 5 X IQi" Hz - . (40) 

We show in Fig. [3] the power spectrum of gravitational 
waves as a function of (comoving) wavenumber k/m. We 
have used different lattices in order to have lattice ar- 
tifacts under control, specially at late times and high 
wavenumbers. We made sure by the choice of lattice 
size and spacing (i.e. fcmin and fcmax) that all relevant 
scales fitted within the simulation. Note, however, that 
the lower bumps are lattice artifacts, due to the physical 
cutoff imposed at the initial condition, that rapidly dis- 
appear with time. We have also checked that the power 
spectrum of the scalar fields follows turbulent scaling af- 
ter mt ~ O(IOO), and we can thus estimate the subse- 
quent evolution of the energy density distributions be- 
yond our simulations. 

IV. LATTICE SIMULATIONS 

The problem of determining the time evolution of a 
quantum field theory is an outstandingly difficult prob- 
lem. In some cases only a few degrees of freedom are 
relevant or else perturbative techniques are applicable. 
However, in our particular case, our interests are focused 
on processes which are necessarily non-linear and non- 
perturbative and involve many degrees of freedom. The 
presence of gravitational fields just contributes with more 
degrees of freedom, but does not complicate matters sig- 
nificantly. 

The lattice formulation allows a first principles ap- 
proach to non-perturbative quantum field theory. The 
existing powerful lattice field theory numerical meth- 
ods rest on the path integral formulation in Euclidean 



space and the existence of a probability measure in 
field space [4l| . However, the problem we are inter- 
ested in is a dynamical process far from equilibium, and 
the corresponding Minkowski path integral formulation 
is neither mathematically well founded nor appropriate 
for numerical studies. There are a series of alterna- 
tive non-perturbative methods which different research 
groups have used to obtain physical results in situa- 
tions similar to ours. These include Hartree's approxi- 
mations [12] togo beyond perturbation theory or large 
N techniques [ij, SI] . It is clear that it is desirable to 
look at this and similar problems with all available tools. 
In the present paper we will use an alternative approx- 
imation to deal with the problem: the classical approx- 
imation. It consists of replacing the quantum evolution 
of the system by its classical evolution, for which there 
are feasible numerical methods available. The quantum 
nature of the problem remains in the stochastic char- 
acter of the initial conditions. This approximation has 
been used with great success by several groups in the 
past [H, [sl] . The advantage of the method is that it is 
fully non-linear and non-perturbative, allows the use of 
gauge fields [3l|, H^l and gives access to the quantities we 
are interested in. 

The validity of the approximation depends on the loss 
of quantum coherence in the evolution of the system. In 
previous papers we studied this problem both in the ab- 
sence of and with gauge fields [2^, HH, [111 . We started 
the evolution of the system at the critical time tc, corre- 
sponding to the end of inflation ig, at which the effective 
mass of the Higgs vanishes, putting all the modes in its 
(free field) ground states. If the couplings are small, since 
the quantum fluctuations of the value of the fields are not 
too large, the non-linear terms in the Hamiltonian of the 
system can be neglected. Then the quantum evolution 
is Gaussian and can be studied exactly. The Hamilto- 
nian has nonetheless a time-dependence coming through 
the time-dependence of the inflaton homogeneous mode. 
This time dependence can always be taken to be linear 
for a sufficiently short time interval after the critical time. 
As a consequence, the dynamics of the eigenmodes dur- 
ing this initial phase differs significantly from mode to 
mode. Most of them have a characteristic harmonic oscil- 
lator behaviour with a frequency depending on the mode 
in question. In the case of the Higgs field, the long-wave 
modes become tachyonic. By looking at expectation val- 
ues of products of these flelds at different times, one re- 
alises that after a while these modes behave and evolve 
like classical modes of an exponentially growing size. The 
process is very fast and therefore the remaining harmonic 
modes can be considered to have remained in their initial 
ground state. 

The fast growth in size of the Higgs field expecta- 
tion value boosts the importance of non-linear terms and 
eventually drives the system into a state where the non- 
linear dynamics, including the back- reaction to the in- 
flaton field, are crucial. For the whole approximation to 
be useful this must happen at a later time than the one 
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in which the low-frequency Higgs modes begin evolving 
as classical fields. In paper I we showed this to be the 
case. Actually, there is a time interval in which classical 
behaviour is already dominant while non-linearities are 
still small. We tested that our results, in the absence 
of gauge fields, were insensitive to the matching time, 
provided it lies within this window. 

The whole idea can then be summarised as follows: 
the tachyonic quantum dynamics of the low momentum 
Higgs modes drive them into classical field behaviour and 
large occupation numbers before the non-linearities and 
back-reaction begin to play a role. It is the subsequent 
non-linear classical behaviour of the field that induces 
the growth of classical inflaton and gravitational field 
components also at low frequencies. It is obvious that 
the quantum nature of the problem is still crucial if one 
studies the behaviour of high momentum modes which 
have low occupation numbers. 

In the present paper we apply the same idea in the 
presence of (gravitational wave) tensor fields. The initial 
quantum evolution of tensor fields is also relatively slow, 
since there are no tachyonic modes. Therefore, it is as- 
sumed not to affect substantially the initial conditions of 
the classical system. 



A. Initial conditions 

Our approach to the dynamics of the system is to as- 
sume that the leading effects under study can be well- 
described by the classical evolution of the system. The 
justification of this point, as explained in the previous 
section, rests upon the fast quantum evolution of the 
long wavelength components of the Higgs field during 
the initial stages after the critical point. All the other 
degrees of freedom will evolve slowly from their initial 
quantum vacuum state. For the Higgs field, the lead- 
ing behaviour is the exponential growth of those modes 
having negative effective mass-squared. The quantum 
evolution of such modes drives the system into a quasi- 
classical regime. It is essential that this regime is reached 
before the non-linearities couple all degrees of freedom to 
each other and questions like back-reaction start to affect 
the results. It is then assumed that it is the essentially 
classical dynamics of that field what matters, and that all 
the long-wavelength components of the inflaton and the 
gauge fields produced by the interaction with the Higgs 
field behave also as classical fields. Of course, this can 
hardly be the case for shorter wavelengths which stay in a 
quantum state with low occupation numbers. However, 
as we can see in Fig. [31 for the range of times studied 
in this paper, the effect of shorter wavelengths is rela- 
tively small, and the spectrum of modes remains always 
dominated by long- wavelengths. 

The full non-linear evolution of the system can then 
be studied using lattice techniques. Our approach is to 
discretize the classical equations of motion of all fields 
in both space and time. The time-like lattice spacing at 



must be smaller than the spatial one ag for the stability of 
the discretized equations. In addition to the ultraviolet 
cut-off one must introduce an infrared cut-off by putting 
the system in a box with periodic boundary conditions. 
We have studied 64^, 128^ and 256^ lattices. Computer 
memory and CPU resources limit us from reaching much 
bigger lattices. Nonetheless, in the spirit of paper I, there 
are a number of checks one can perform to ensure that our 
results are physical and are not biased, within errors, by 
the approximations introduced, see Fig. [3] Our problem 
has several physical scales which control different time- 
regimes and observables. Thus, it is not always an easy 
matter to place these scales in the window defined by 
our ultraviolet and infrared cut-offs. For example, in 
addition to the Higgs and inflaton mass there is a scale M 
associated to the inflaton velocity which is particularly 
relevant in determining the bubble sizes and collisions. 
Then, when we want to study a stage of the evolution 
in particular, we make the selection of the volume and 
cutoff most suitable. 

Since in this paper we are more interested in under- 
standing the phenomenon of GW production, rather than 
concentrating in a particular model, our attitude has 
been to modify the parameters of the model in order to 
sit in a region where our results are insensitive to the cut- 
offs. This is no doubt a necessary first step to determine 
the requirements and viability of the study of any partic- 
ular model. In particular, we have thouroughly studied 
a model with — 2X — 1/4, but we have checked that 
other values of the parameters do not change our results 
significantly. 

The initial conditions of the fields follow the prescrip- 
tion from paper I. The Higgs modes (f>k are solutions of 
the coupled evolution equations, which can be rewrit- 
ten as (j)'^ + {k^ — T)(j)k = 0, with r = M{t — tc) and 
M = {2Vy/'^m. The time-dependent Higgs mass follows 
from the initial inflaton field homogeneous component, 
XoiU) = Xc(l - Vm{U - tc)) and Xq{U) = ~XcVm. The 
Higgs modes with k/M > ^Jri are set to zero, while the 
rest are determined by a Gaussian random field of zero 
mean distributed according to the Rayleigh distribution 

P{\cj)k\)d\(l>k\dek = exp — — — , (41) 

with a uniform random phase 9^ £ [0, 27r] and disper- 
sion given by cr^ = = P{k,Ti)/k^, where Pik.n) 
is the power spectrum of the initial Higgs quantum fiuc- 
tuations in the background of the homogeneous infiaton, 
computed in the linear approximation. In the region of 
low momentum modes it is well approximated by 

2fcPapp(fc,rO-fc3(l + yl(r,)fc'e-^(^')"') , (42) 

where A{Ti) and B{Ti) are parameters extracted from a fit 
of this form to the exact power spectrum given in paper I. 
In the classical limit, the conjugate momentum 4>k{T) 
is uniquely determined through <f)k{T) ~ F(k,T)(f>k{T), 
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FIG. 4: The tachyonic growth of the Higgs' spectrum, from 
mt = 5 to mt — 10. We compare simulations of different 
sizes (pmin = 0.01 — 0.03) and A'^ = 256, with the anayhtical 
expressions. 

where F{k,T) = Im(i/fe(T)/*(T))/|/fc(r)|2, see paper I. 
In the region of low momenta, F(k, Ti) can be well aprox- 
imated by 



F{k,n) = 



[l + A(T,)e-^(-')'=']' 



(43) 



where A{Ti) and B{Ti) are the previous coefficients for the 
amplitude of the field fluctuations, while C{Ti) and D{Ti) 
are new coefficients obtained fitting the exact expression 

oiF{k,n). 

The rest of the fields (the inflaton non-zero modes and 
the gravitational waves), are supposed to start from the 
vacuum, and therefore they are semiclassically set to zero 
initially in the simulations. Their coupling to the Higgs 
modes will drive their evolution, giving rise to a rapid (ex- 
ponential) growth of the GW and inflaton modes. Their 
subsequent non-linear evolution will be well described by 
the lattice simulations. 

In the next subsections we will describe the different 
evolution stages found in our simulations. 



B. Tachyonic growth 

In this subsection we will compare the analytical es- 
timates with our numerical simulations for the initial 
tachyonic growth of the Higgs modes and the subsequent 
growth of gravitational waves. The first check is that the 
Higgs modes grow according to paper I. There we found 
that 



with A[t) and B{t) are given in paper I, 



(44) 



2/1 /Q^2/3 

A[r) = -1^^^ Bi2(r) , B{t) ^ - 1) , (45) 

which are valid for r > 1, and where Bi(z) is the Airy 
hmction of the second kind. Indeed, we can see in Fig. 3] 
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FIG. 5: The Fourier transform of the anisotropic stress tensor. 
We compare the numerical simulations of nii(fc, t) for Pmin = 
0.01 with the analytical expressions (dashed lines) for mt = 
5— 10, i.e. during the tachyonic growth. The small deviations 
at fc < m are simulation artifacts due to the initial UV cut-off 
implementation and soon disappear. 

that the initial growth, from mt = 6 to mt = 10, fol- 
lows precisely the analytical expression, once taken into 
account that in Eq. (|44p the wavenumber k and time r 
are given in units of M = {2VY^^m. 

The comparison between the tensor modes hij{k,t) 
and the numerical results is somewhat more complicated. 
We should ffi'st compute the effective anisotropic ten- 
sor Tij (k, t) ([M)) from the gradients of the Higgs field 
(those of the inflaton are not relevant during the tachy- 
onic growth), as follows. 



H.,(k,i) = 
where 



(27r)3/2 



(27r) 



3/2 



2qx 



(46) 



(47) 



After performing the integral in x and using the delta 
function to eliminate q', we make a change of variables 
q q-f k/2, and integrate over q, with which the Fourier 
transform of the anisotropic stress tensor becomes 



H/j" (k, t) — ki kj 



A{r) 
B{t)V2 



1 B(t)P 
- 0- — 

2' ' 4 



(48) 



which gives a very good approximation to the numerical 
results, see Fig. El with ^'(1/2,0; z) ~ (tt^^ +2)^^/2 being 
the Kummer function. 

Finally, with the use of Hij (k, t) , we can compute the 
tensor fields. 



sink{t~t') ~ 



(49) 



h,j{k,t) = (167rG) / dt' 
Jo 

dQh,j{k,t) = (167rG) / dt' cos k{t - t')Ilij. (50) 
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Using the analytic expression in Eq. (|48p one can perform 
the integrals and obtain expressions that agree surpris- 
ingly well with the numerical estimates. This allows one 
to compute the density in gravitational waves, Pg„, at 
least during the initial tachyonic stage in terms of ana- 
lytical functions, and we reproduce the numerical results, 
see Fig. H 

We will now compare our numerical results with the 
analytical estimates. The tachyonic growth is dom- 
inated by the faster than exponential growth of the 
Higgs modes towards the true vacuum. The (traceless) 
anisotropic strees tensor 11^ grows rapidly to a value of 
order fc^|0p ^ lO^^m^u^, which gives a tensor pertur- 
bation 



1/2 



167rGu^(mAi)2lO- 



and an energy density in GW, 



(51) 



(52) 



for niAt ~ 16. In the case at hand, with v = 10^^ Mp, 
we find p^^/po ~ 10~^ at symmetry breaking, which 
coincides with the numerical simulations at that time, 
see Fig. H 

As shown in Ref. [1^, the spinodal instabilities grow 
following the statistics of a gaussian random field, and 
therefore one can use the formalism of |45|] to estimate the 
number of peaks or lumps in the Higgs spatial distribu- 
tion just before symmetry breaking. As we will discuss in 
the next section, these lumps will give rise via non-linear 
growth to lump invagination and the formation of bubble- 
like structures with large density gradients, expanding at 
the speed of light and colliding among themselves giving 
rise to a large GWB. The size of the bubbles upon col- 
lision is essentially determined by the distance between 
peaks at the time of symmetry breaking, but this can be 
computed directly from the analysis of gaussian random 
fields, as performed in Ref. 29]. 

This analysis works only for the initial (linear) stage 
before symmetry breaking. Nevertheless, we expect the 
results to extrapolate to later times since once a bubble is 
formed around a peak, it remains there at a fixed distance 
from other bubbles. This will give us an idea of the size 
of the bubbles at the time of collision. 

We estimate the number density of peaks as [45| 



'^poak(T) 



,^o(T)-3(i''-l)exp[-z.V2], (53) 



where v = 4>c/<j{t) corresponds to the ratio of the field 
threshold 0c over the dispersion 



(54) 



with A{t) and B{t) given in Eq. (^5)) . The average size 
of the gaussian lumps is £,o{t) — 2i?^/^(T) to~\ where 
time is given in units r = (2V)^/'^ mt, see Ref. [29| 



The distance between peaks can be estimated as twice 
the radius of the average bubble, with volume Vpoak = 
47r/3 i?p^^j^. Since the total volume is divided into 
A'peak bubbles, we find 



^pcak 



2R 



peak 



^( 

ma V 



6 



1/3 



(55) 



which is typically of order 30 to 40 lattice units, for (pc — 
0.5 - 0.8, V = 0.024 and A = 0.125, with lattices sizes 
given by Pmin = 0.15 m and N = 128. 

What is interesting is that decreasing either A or V, the 
distance between initial lumps increases and thus also the 
size of the final bubbles upon collision. As we will show 
in the next section, the amplitude of GW depends on the 
bubble size squared, and therefore it is expected that for 
lower lambda we should have larger GW amplitude. We 
have not seen, however, such an increase in amplitude, 
but a detailed analysis is underway and will be presented 
elsewhere. 



C. Bubble collisions 

The production of gravitational waves in the next stage 
proceeds through 'bubble' collisions. In Ref. [l^ we 
showed explicitly that symmetry breaking is not at all 
a homogeneous process. During the breaking of the 
symmetry, the Higgs field develops lumps whose peaks 
grow up to a maximum value |0|max/v = 4/3 (see pa- 
per I), and then decrease creating approximately spheri- 
cally symmetric bubbles, with ridges that remain above 
1 01 — V. Finally, neighboring bubbles coUide and the 
symmetry gets further broken through the generation of 
higher momentum modes. Since initially only the Higgs 
field sources the anisotropic stress-tensor Tlij, then we 
expect the formation of structures (see section IV. A) in 
the tensor metric perturbation, correlated with the Higgs 
lumps. The dependence of the hij tensor on the gradient 
of the Higgs field, see Eq. p^ . is responsible of the for- 
mation of those structures in the energy density spatial 
distribution of the GWB. 

In section V of this paper we will give account of the 
explicit form of the structures developed in the spatial 
distribution of related with the first collisions among 
the bubble-like structures of the Higgs field. We will 
present simultaneously the evolution of both the Higgs' 
spatial distribution when the first bubbles start colliding, 
and of the corresponding structures in the GW energy 
density p^^ . We leave for a forthcoming publication the 
details of an analytical formalism describing the forma- 
tion and subsequent evolution of such GW structures. 
In this sub-section we will just give an estimate of the 
burst in GW produced by the first collisions of the Higgs 
bubble-like structures. 

Thus, as described in Ref. [l^l for the case of first order 
phase transitions, which involves the collision of vacuum 
bubbles, we can give a simple estimate of the order of 
magnitude of the energy fraction radiated in the form 
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of gravitational waves when two Higgs bubble-like struc- 
tures collide. In general, the problem of two colliding 
bubbles has several time and length scales: the duration 
of the collision. At; the bubbles' radius R at the moment 
of the collision; and the relative speed of the bubble walls. 
In section IV. B we found that the typical size of bubbles 
upon collisions, is of the order of i? « 10m~^, while the 
growth of the bubble's wall is relativistic. Then we can 
assume than the time scale associated with bubble colli- 
sions is also At ~ R. Assuming the bubble walls contain 
most of the energy density, and since they travel close to 
the speed of light, see paper I, it is expected that the 
asymmetric collisions will copiously produce GW. 

Far from a source that produces gravitational radia- 
tion, the dominat contribution to the amplitude of GW 
is given by the acceleration of the quadrupole moment 
of the Higgs field distribution. Given the energy den- 
sity of the Higgs field, pfji we can compute the (re- 
duced) quadrupole moment of the Higgs field spatial 
distribution, Qij — J (Px(xiXj — x^Sij /3) pB_ix), such 
that the amplitude of the gravitational radiation, in the 
TT gauge, is given by hij ~ {2G/r)Qij. A signifi- 
cant amount of energy can be emitted in the form of 
gravitational radiation whenever the quadrupole moment 
changes significantly fast: through the bubble collisions 
in this case. The power carried by these waves can be 
obtained via (|29p as 

^ow = ^ / dn (q./Q'') . (56) 

Omitting indices for simplicity, as the power emitted in 
gravitational waves in the quadrupole approximation is 
of order ~ G{Q)^, while the quadrupole moment is 
of order Q ~ R^pa, we can estimate the power emitted 
in GW upon the collision of two Higgs bubbles as 

^ow ^ GplR' (57) 

The fraction of energy density carried by these waves, 
Pgw ~ Pcw^t/R^ ~ Pg^/R^ ~ Gp^R^, compared 
to that of the initial energy stored in the two bubble-like 
structures of the Higgs field, will be p^^/pn = GpuR^- 
Since the expansion of the universe is negligible during 
the bubble collision stage, the energy that drives infla- 
ton, po ^ m^ti^, is transferred essentially to the Higgs 
modes during preheating, within an order of magnitude, 
see Fig. [2] Thus, recalling that R ^ 10m~^, the total 
fraction of energy in GW produced during the bubble 
collisions to that stored in the Higgs lumps formed at 
symmetry breaking, is given by 

^ r^O.lGpoR^ {v/Mp)\ (58) 
Pa 

giving an amplitude which is of the same order as is ob- 
served in the numerical simulations, see Fig.[2l Of course, 
an exhaustive analytical treatment of the production of 
GW during this stage of bubble collisions remains to be 
done, but we leave it for a future publication. 
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FIG. 6: Variance of the Inflaton and the Higgs field as a 
function of time, the former normalized to its critical value, 
the latter normafized to its v.e.v.. As expected in a turbulent 
regime, these variances follow a power law t~'^^ with p a 
certain critical exponent, although the slope of the Inflaton's 
variances evolves in time. The curves are produced from an 
average over 10 different statistical realizations. 

D. Turbulence 

The development of a turbulent stage is expected from 
the point of view of classical fields, as turbulence usu- 
ally appears whenever there exists an active (stationary) 
source of energy localized at some scale in Fourier 
space. As first pointed out by Ref. in reheating sce- 
narios the coherently oscillating inflaton zero-mode plays 
the role of the pumping-energy source, acting at a well 
deflned scale fcjn in Fourier space, given by the frequency 
of the inflaton oscillations. Thus, the inflaton zero-mode 
pumps energy into the rest of the fields that couple to it 
as well as into the non-zero modes of the inflaton field 
itself. Apart from kin, there is no other scale in Fourier 
space where energy is accummulated, dissipated and/or 
infused. So, as turbulence is characterized by the trans- 
port of some conserved quantity, energy in our case, we 
should expect a flow of energy from k-m towards higher 
(direct cascade) or smaller (inverse cascade) momentum 
modes. In typical turbulent regimes of classical fluids, 
there exits a sink in Fourier space, corresponding to that 
scale at which the (direct) cascade stops and energy gets 
dissipated. However, in our problem there is no such sink 
so that the transported energy cannot be dissipated, but 
instead it is used to populate high-momentum modes. 
For the problem at hand, there exists a natural initial 
cut-off fcout ~ X^^'^v, such that only long wave modes 
within k < fcout, develop the spinodal instability. Even- 
tually, after the tachyonic growth has ended and the first 
Higgs' bubble-like structures have collided, the turbulent 
regime is established. Then the energy flows from small 
to greater scales in Fourier space, which translates into 
the increase of kout in time. 

In Ref. [3^ we already accounted for the turbulent 
stage reached in a hybrid model with gauge fields. How- 
ever, we don't consider gauge fields here, so the number 
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FIG. 7: Some snapshots of the evolution of the spectral par- 
ticle occupation numbers of the Higgs field at different times, 
each averaged over 10 statistical realizations. We multiply 
them by fc* so we can see better the scaling behaviour. In 
the upper right corner, we plot the inverse relation of (I60|l . 
no{kt~^) = t~'^n{k,t), also averaged over 10 realizations for 
each time. The scaling behaviour predicted by wave kinetic 
turbulent theory '461], is clearly verified. 
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FIG. 8: Different times of the evolution of the particle occu- 
pation numbers spectra of the Infiaton, multiplied by fc*, and 
averaged over 10 statistical realizations or each time. Again, 
in the upper right corner, we plot the inverse relation of (|60[1 . 
no(fc) = t''^n{kt^ ,t), also averaged over 10 realizations for 
each time. 

of degrees of freedom is different from that of Ref. (s^l 
and, therefore, the turbulent dynamics of the Infiaton 
and the Higgs fields should be different. In particular, 
when the turbulence has been fully established, if the 
wave (kinetic) turbulence regime of the fields' dynamics 
is valid, the time evolution of the variance of a turbulent 
field /(x, t), should follow a power-law-like scaling [4^ 

Var(/(<)) = - {f[t)f (X t-^'P , (59) 

withp = l/(2iV— 1) and N the number of scattering fields 
in a 'point-like collision'. In fact, such time behaviour 
corresponds only to the case of the so called free turbu- 
lence, when the energy stored in the pumping source is 
subdominant to the energy in the turbulent fields. In our 
case, this condition is reached very soon after the sym- 
metry breaking, so we don't expect a significant stage 



of driven turbulence, which would make the variance to 
increase (Only the infiaton seems to increase its variance 
between mt — 10 and mt — 30, but it is not very pro- 
nounced) . In Fig. [HI we have plotted the time evolution 
of the variances of the Infiaton x and of the Higgs mod- 
ulus (j) ~ ■\/Sa"Wi a'^^ fitted the data with a power-law 
like ([59|l . obtaining 

Infiaton: p-^ = 5.1 ± 0.2, [35:85] 

Infiaton: p-^ = 9.03 ± 0.03, [350:2000] 

Higgs; = 7.02 ± 0.01, [50:2000] 

where the last brackets on the right correspond to the 
range in time (in units of m~^) for which we fitted the 
data. As can be seen in Fig. [SI the slope of the Higgs 
field (in logarithmic scale), 2pjj ^ 2/7, remains approxi- 
mately constant in time, corresponding to a 4-field dom- 
inant interaction. However, the slope of the Infiaton's 
variance increases in time, i.e. the critical exponent p^ of 
the Infiaton decreases, until it reaches a stationary stage 
at mt ^ 100. Since Pj is related to the number N of fields 
interacting in a collision, if there was a change from one 
dominant multi-field interaction to another, this should 
produce a time-dependent effective p^ , as seen in Fig. [G] 
However, we will not try to explain here the origin of 
such an effective critical exponents as extracted from the 
simulations. We will just stress that we have checked 
the robustness of those values under different lattice con- 
figurations (A^, Pmin) and different statistical realizations, 
discarding this way a possible lattice artefact effect. As 
we will see, the critical exponents p determines the speed 
with which the turbulent particle distribution moves over 
momentum space, so this is a crucial parameter. More- 
over, although both the classical modes of the Infiaton 
and the Higgs contribute to the production of GW, the 
Infiaton's occupation numbers decrease faster than those 
of the Higgs so, after a given time, only the Higgs' modes 
remain as the main source of GW. 

Actually, when turbulence is developed, it is expected 
that the distribution function of the classical turbulent 
fields, the infiaton and the Higgs here, follow a self-similar 
evolution 

n{k,t) =t-'^PnQ{kt-P) (60) 

with p the critical exponent of the fields' variances and 
7 a certain factor ~ 0(1), which depends on the type of 
turbulence developed. It is precisely this way that the 
exponent p determines the speed of the particles' distri- 
bution in momentum space: given a specific scale kc such 
that, for example, the occupation number has a maxi- 
mum, that scale evolves in time as kc{t) = kc{t()){t/to)P . 
We have seen that the evolution of the Higgs occupation 
number follows Eq. (|60p with p « 1/7, as expected from 
the Higgs variance, and 7 « 2.7. Whereas the evolution 
of the Infiaton occupation number follows ((60)) even more 
accurately than the Higgs, with an "effective" exponent 
p « 1/5, and 7 w 3.9. Since the slope of the infiaton's 
variance changes in time, the value of the exponents of 



15 



1e-05 




FIG. 9: Time evolution of the GW spectra from mt = 6 to 
mt = 2000. The amplitude of the spectra seems to saturate 
after mt ~ 100, although the high momentum tail still moves 
slowly to higher values of k during the turbulent stage. 

the inflaton's scaling relation will require further investi- 
gation. However, despite this time evolution of the Infla- 
ton variance, Eq. ([TO]) is very well fulfilled by the Inflaton 
with the given effective exponents. So we can perfectly 
obtain the universal rio(fc) function for the Inflaton as 
well as for the Higgs. 

In Figs. [7] and [5] we have plotted the occupation num- 
bers of the Higgs and the Inflaton, also inverting the re- 
lation of Eq. (j60p in order to extract the universal time- 
independent no(fc) functions of each field. As shown in 
those figures, the distributions follow the expected scal- 
ing behaviour. However, for the range of interest of fc, 
there are small discrepancies of order 0.1-4% (depending 

on k) among the different universal functions riQ^k), as 
obtained inverting Eq. (|60p at different times mti. The 
universal functions n(){k) plotted in Figs. [7] and [8] have 
been obtained from averaging over ten statistical realiza- 
tions for each time. 

The advantage of the development of a turbulence be- 
haviour is obvious: it allows us to extrapolate the time 
evolution of the fields' distributions till later times be- 
yond the one we can reach with the simulations. More- 
over, the fact that the turbulence develops so early after 
the tachyonic instability, also allow us to check for a long 
time of the simulation, the goodness of the description 
of the dynamics of the fields, given by the turbulent ki- 
netic theory developed in Ref. ,4^. We have fitted the 
averaged universal functions no{k) with expressions of 
the form k^ no{k) = P{k)e-Q'^''\ with P{k) and Q{k) 
polynomials in fc, giving: 

Inflaton : P{k) = 486.2fc3 Q{k) = 6.39fc 

(61) 

Higgs : Pik) = 2.96fc3 Q{k) = 2.26P - 3.18k 

There is no fundamental meaning for these expressions, 
but it is very useful to have such an analytical control 
over no(fc), since this allows us to track the time-evolution 



of n{k, t) through Eq. (pD)) . Actually, the classical regime 
of the evolution of some bosonic fields ends when the 
system can be relaxed to the Bose-Einstein distribution. 
We are now going to estimate the moment in which the 
initial energy density gets fully transferred to the Higgs 
classical modes. Using Eg. ([60)) and the fit ([6T|) . we find 
that the initial energy density is totally transfered to the 
Higgs when (in units m = 1) 

1 f dk k"^ , , T.565 /^^\ 
Po^:^=y^^fcn(fc,^) = ^.(-^)^ (62) 

where we have assumed that the Higgs' modes have en- 
ergy Ek{k,t) = kn{k,t). In our case, with A — 1/8, 
the conversion of the initial energy density into Higgs 
particles and therefore into radiation is complete by 
mt ~ 6 X 10'*. Therefore, if we consider this value as 
a lower bound for the time that classical turbulence re- 
quires to end, we see that turbulence last for a very 
long time compared to the time-scale of the initial tachy- 
onic and bubbly stages. Thus, if GW were significatively 
sourced during turbulence, one should take into account 
corrections from the expansion of the universe. 

In Fig. [HI we show the evolution of the GW spectra up 
to times mt = 2000, for a lattice of (N,pniin) = (128, .15). 
It is clear from that figure that the amplitude of the GW 
saturates to a value of order //9o « 2 • 10"^ At mt 
~ 50, the maximum amplitude of the spectra has al- 
ready reached Pgw/po ~ 10~^, while at time mt ~ 100, 
the maximum has only grown a factor of 2 with respect 
to mt « 50. From times mt « 150 till the maximum 
time we reached in the simulations, mt = 2000, the max- 
imum of the amplitude of the spectrum does not seem to 
change significantly, slowly increasing from « 2 • 10 ~^ to 
« 2.5 • 10^^. Despite this saturation, we see in the simu- 
lations that the the long momentum tail of the spectrum 
keeps moving towards greater values. This displacement 
is precisely what one would expect from turbulence, al- 
though it is clear that the amplitude of the new high mo- 
mentum modes never exceed that of lower momentum. 
In order to disscard that this displacement towards the 
UV is not a numerical artefact, one should further inves- 
tigate the role played by the turbulent scalar fields as a 
source of GW. Here, we just want to remark that the tur- 
bulent motions of the scalar fields, seem not to increase 
significatively anymore the total amplitude of the GW 
spectrum. Indeed, in a recent paper [25| where GW pro- 
duction at reheating is also considered, it is stated that 
GW production from turbulent motion of classical scalar 
fields, should be very supressed. That is apparently what 
we observe in our simulations although, as pointed above, 
this issue should be investigated in a more detailed way. 
Anyway, here we can conclude that the expansion of the 
Universe during reheating in these hybrid models, does 
not play an important role during the time of GW pro- 
duction, and therefore we can be safely ignore it. 
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FIG. 10: Model: A = 10"^ g"^ = 1. In the left picture we show a spatial section of (0). We can see how a spherical lump is 
growhing. In the right picture we can see the structure of the paw in the same place. A ring is forming around the Higgs lump. 
More complex structures are formed in the regions in which the Higgs bubbles are next, and the GW grow in the boundary of 
this lumps, where the gradient of the Higgs and therefore the Ilij tensor grow in this region. 



V. SPATIAL SECTIONS AND LOCAL GW 
PRODUCTION 

In this section, we show a sequence of snapshots {mt — 
5 — 20) of the evolution of the spatial distribution, before 
the fields are driven to the turbulent stage. We find that 
the first stages of the GW dynamics is strongly correlated 
with the dynamics of the Higgs oscillations that give rise 
to symmetry breaking. A qualitative way of understand- 
ing this question is to analyse the spatial structure of 
the liij tensor, built from spatial gradients of the Higgs 
and inflaton fields. Since the oscillations of (</>) are due to 
rapid changes of the Higgs' values in its way of symmetry 
breaking, this induces great variations in the behaviour 
of the spatial gradients. We are now going to analyse 
briefiy the different stages showing the most representa- 
tives images. Shortly, it will be available in our web page, 
a bigger selection of pictures and movies [47 ] . 

An interesting conclusion from the set of Figs.fTUl — [HI 
is that the Higgs evolution from lump growth, through 
invagination to bubble collisions, has a direct transla- 
tion into the corresponding growth of gravitational wave 
energy density. Not only does the volume-averaged am- 
plitude pgw follow the Higgs time evolution, but the in- 
dividual local features in the GWB seem to correspond 
very closely with the Higgs features. 

In the first stage both Higgs and GW backgrounds 
grow very fast. The lumps which grow in the Higgs back- 
ground induce structures around these, through the gra- 
dients appearing in the H^j tensor. The geometry of the 
gravitational structures comes from the position of the 
Higgs lump. A typical structure for an isolated lump is 



a ring of gravitational waves, see Fig. [TOl More complex 
structures can be formed around the position of the Higgs 
lumps. Before symmetry breaking these lumps grow ac- 
cording to the previous analysis, generating domains with 
a great density of gravitational energy. 

The second stage begins when (0) = v and the sym- 
metry breaking starts, then the Higgs lumps invaginate 
and expand, producing the growth of gravitational waves 
around of these structures, see Figs. fTT]and[T2l one can 
see that whenever the bubble walls expand, the variation 
in the gradient of the Higgs' field induces the expansion 
of the GW ring. In the case of the rings, if the lump 
is very isolated, the expansion induces the ring to dilute 
and disappear, by Gauss law. In practice, however, the 
lumps are never isolated and bubbles collide before the 
gradients (and thus the GW) die away. 

In the case when two Higgs' bubble-like structures are 
close by, the expansion of their walls compresses the GW 
structures. This expansion continues until the first Higgs 
oscillation, see Fig. [3l If the distance between Higgs' 
structures is small, then the GW can be diluted, whereas 
in the other remnant string-like GW structure sur- 

vives, and when the Higgs background goes to zero this 
GW structure becomes divided into two waves that prop- 
agate in opposite directions, as one can see in Figs. [T2] 
and [Ml which show four snapshots of this process. A 
similar behaviour is observed in the second oscillations. 

Finally, the wave fronts propagate, colliding among 
themselves, driving the system to the stage of turbulence. 
We will leave for a future publication the detailed anal- 
ysis of the GW production at the bubble collisions and 
the subsequent turbulent period. 
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FIG. 11: Here we have got the time evolution of the previous 
ring (picture nop near the symmetry breaking. The bubble is 
growing (mt = 10-11), until the symmetry breaking time (mt 
= 12). 



FIG. 12: The Higgs lump begins to invaginate, and the GW 
ring expands (mt = 12-13). A similar behavior is observed 
in a smaller lump below the biggest Higgs lump, in the same 
pictures. 



VI. GRAVITATIONAL WAVES FROM 
CHAOTIC INFLATION 

The production of a relic GWB at reheating was 
first addressed by Khlebnikov and Tkachev (KT) in 
Ref. [27\, both for the quadratic and quartic chaotic in- 
flation scenarios. In these models, the long-wavelength 
part of the spectrum is dominated by the gravitational 
bremsstrahlung associated with the scattering of the ex- 
tra scalar particles off the inflaton condensate, 'evaporat- 
ing' this way the inflaton particles. Using this fact, KT 
estimated analytically the amplitude of the power spec- 
tra of GW for the low frequency end of the spectrum, 
corresponding to wavelengths of order the size of the hori- 
zon at rescattering. Moreover, KT also studied the GW 
power spectra numerically, although just for the massless 
inflaton case. Recently, chaotic scenarios were revisited 
in Ref. [23, 24], accompanied by more precise numerical 
simulations at different energy scales, including the case 
of a massive inflaton. Finally, also very recently, Ref. [2^ 
studied in a very detailed way, both analitical and numer- 
ically, the evolution of GW produced at preheating in the 
case of a massless inflaton with an extra scalar field. 

In Refs. [llj and [11], the procedure to compute the 
GW from reheating relied on Weinberg's formula for the 
energy carried by a weak gravitational radiative field in 
flat space-time [4^. However, in chaotic models, the 



expansion of the universe can not be neglected during 
reheating, so Weinberg's formula can only be used in 
an approximated way, if the evolution of the universe 
is considered as an adiabatic sequence of stationary uni- 
verses. Rescaling fields by a conformal transformation, 
their evolution equations can be solved with a numerical 
integrator, while the evolution of the scale factor can be 
calculated analytically. Discretizing the time, the phys- 
ical variables can be recovered from the conformal ones 
in each time step, thus allowing to compute the energy 
of gravitational waves in terms of the physical fields. In 
this paper, however, we adopt another approach^ that 
takes into account expansion of the universe in a self- 
consistent manner, and let us calculate at any time the 
energy density and power spectra of the GW produced 
at reheating. As explained in section III and applied to 
the case of hybrid inflation in sections IV and V, we just 
solve numerically Eq. together with those eqs. of 

the other Bose fields and the scale factor, Eqs. (iQjl. pO]) 
and Eqs.(ini),([ni). Then, using the projector (HH) into 



^ Note that Refs. [24l | and [SH ] also work in the same theoretical 
framework, considering the TT tensor perturbations on top of 
a flat FRW space. However, we use a different way to extract 
numerically the GW power spectra, relying on the conmutating 
procedure, as detailed explained in subsection III.B 
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FIG. 13: After symmetry breaking the expansion of Higgs 
lump compresses the GW, until the Higgs gradient changes 
in the first oscillation (mt =13-13.5). 



FIG. 14: At this moment, when the Higgs falls, the GW 
strutrure is divided in two waves (mt = 14). These wave 
fronts propagates in opposite directions (mt — 14.5). 



the (Fourier transformed) solution of Eq. (|23p , we recover 
the TT d.o.f corresponding to GW. This way, we can 
monitor the total energy density in GW using Eq. 
or track the evolution of the power spectrum. Using this 
technique, we will show in this section that we repro- 
duce, for specific chaotic models, similar results to those 
of other authors. 

We adapted the publicly available LATTICEEASY 
code [1^ , taking advantage of the structure of the code 
itself, incorparating the evolution of Eq. p4)) . together 
with the equations of the scalar fields, Eqs. (O and PU)) . 
into the staggered leapfrog integrator routine. This way, 
we can solve at the same time the dynamics of the scalar 
and tensor fields, within the framework of an expanding 
FRW universe Eqs.dH]) and p^ . 

In particular, we will concentrate only in an scenario 
with a massless inflaton x, either accompanied or not by 
an extra scalar field (f). In the following, we will describe 
the numerical results for GW production at reheating in 
such scenarios, described by the potential 



and the physical fields by a conformal transformation as 



Rescaling the time by 



a(0) 



(63) 



(64) 
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(65) 
(66) 



then the equations of motion of the inflaton and of the 
extra scalar field, Eqs. ^ and PT7|) . can be rewritten in 
terms of the conformal variables as 



a 
a 
a" 



Xc - V'xc Xc + ixl + q^Dxc = 



-Xc + qxc(l>c = , 



(67) 
(68) 



where the prime denotes derivative with respect to con- 
formal time. Since the universe expands as radiation-like 
in these scenarios, a(r) ~ t, so the terms proportional 
to a" /a in Eqs. ((67)) and ((68)) are soon zero, as explicitly 
checked in the simulations. Thanks to this, the model is 
conformal to Minkowski. 

The parameter q = /\ controls the strength and 
width of the resonance. For the case of a massless in- 
flaton without an extra scalar field, we just set q = Q 
in Eq. ((^ and ignore Eq. ((S5)) . However, in that case, 
fluctuations of the inflaton also grow via parametric res- 
onance. Actually, they grow as if they were fluctuations 
of a scalar field coupled to the zero-mode of the inflaton 
with effective couplig q = /X = 3, see Ref. f49| . 
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FIG. 15: The spectrum of the gravitational waves' energy den- 
sity, for coupled case with A = 10"" and jX = 120. The 
spectrum is shown accumulated up to different times during 
GW production, so one can see its evolution. At each time, 
it is normalized to the total instant density. This plot corre- 
sponds to a N = 128 lattice simulation, from r = to r = 240. 
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FIG. 16: The spectrum of the gravitational waves' energy 
density, for the pure case, with A = 10"^*. Again, we show 
the spectrum accumulated up to different times during GW 
production, normalized to the total instant density at each 
time. The plot corresponds to a N = 128 lattice simulation, 
from r = to r = 2000. 



Following Refs. [2l| and [H, we set A = 10"" and 
q — 120. Since, this case is also computed in [2^, we 
can also compare our results with theirs. Moreover, we 
also present results for the pure model with no extra 
scalar field, a case only shown in Ref. [^ij. 

We begin our simulations at the end of inflation, when 
the homogeneous inflaton verifies xo ~ 0.342Aip and 
Xo ~ 0. We took initial quantum (conformal) fluctua- 
tions l/\/2fc for all the modes up to a certain cut-off, 
and only added an initial zero-mode for the inflaton, 
Xc(0) = 1, Xc(0)' = 0. In Figs. [I5] and [M we show 
the evolution of f^^,,- during reheating, normalized to 
the instant density at each time step, for the coupled 
and the pure case, respectively. In the case with an extra 
scalar field, the amplitude of the GWB saturates at the 
end of parametric resonance, when the fields variances 
have been stabilized. This is the beginnig of the turbu- 
lent stage in the scalar fields, which seems not to source 
anymore the production of GWs, as already stated in 
Refs.f23l.l2H. For the pure case, we also see the saturation 
of the amplitude of the spectra, see Fig. [111 although the 
long momenta tail seems to slightly move toward higher 
values. 

Of course, in either case, with and without an extra 
field (/), in order to predict today's spectral window of the 
GW spectrum, we have, first, to normalize their energy 
density at the end of GW production to the total energy 
density at that moment. Secondly, we have to redshift 
the GW spectra from that moment of reheating, taking 
into account that the rate of expansion have changed 
significantly since the end of infiation, see Eg. ([55)1 . In 
particular, the shape and amplitude of GW spectra for 
the case with the extra scalar field coupled to the infia- 
ton with q — 120, see Fi g. [ 171 seems to coincide with the 
espectra shown in Ref. [25] . On the other hand, we also 
reproduce in Fig. [T71 a similar spectra to the one shown 
in [2l[, for the case of the pure quartic model. Thanks 
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FIG. 17: Today's ratio of gravitational waves normalized to 
radiation energy density, for both the coupled and the pure 
case. We took Qt/go — 100 to redshift the spectra from the 
time of the end of production till today. 

to the tremendous gain in computer power, we were able 
to resolve the 'spiky' pattern of the spectrum with great 
resolution. For the first time, it is clearly observed the 
exponential tail for large frequencies, see Figs. [THl [T71 
not shown in Ref. 21]. The most remarkable fact, is 
that we also confirm that the peak structure in the GW 
power spectrum, see Fig. I16[ remains clearly visible at 
times much later than the one at which those peaks have 
dissapeared in the scalar fields' power spectrum. So, as 
pointed out in Ref. [2l[ , this characteristic feature distin- 
guish this particular model from any other. 

Let us emphasize that we have run the simulations 
till times much greater than that of the end of the res- 
onance stage, both for the pure and the coupled case. 
The role of the turbulence period after preheating seems, 
therefore, not to be very important, despite its long du- 
ration. Apparently, the no-go theorem about the suppre- 
sion of GW at turbulence, discussed in [2^, is fulfilled. 
In Refs. [31I, [m it was pointed out that gauge couplings 
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or trilinear interactions could be responsible for a fast the resonace stage, in principle this should not affect the 
thermalization of the universe after inflation (see also results shown above. 
Ref. ^52]), but as long as this takes place after the end of 
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FIG. 18: The sensitivity of the different gravitational wave experiments, present and future, compared with the possible 
stochastic backgrounds; we include the White Dwarf Binaries (WDB) ^5(?| and chaotic preheating coupled and pure) for 

comparison. Note the two well differentiated backgrounds from high-scale and low-scale hybrid inflation. The bound marked 
(?) is estimated from ultra high frequency laser interferometers' expectations p^ . 



VII. CONCLUSIONS 

To summarize, we have shown that hybrid models 
are very efficient generators of gravitational waves at 
preheating, in three well defined stages, first via the 
tachyonic growth of Higgs modes, whose gradients act 
as sources of gravity waves; then via the collisions 
of highly relativistic bubble-like structures with large 
amounts of energy density, and finally via the turbulent 
regime (although this effect does not seem to be very 
significant in the presence of scalar sources), which 
drives the system towards thermalization. These waves 
remain decoupled since the moment of their production, 
and thus the predicted amplitude and shape of the 
gravitational wave spectrum today can be used as a 
probe of the reheating period in the very early universe. 
The characteristic spectrum can be used to distinguish 
between this stochastic background and others, like 
those arising from NS-NS and BH-BH coalescence, 
which are decreasing with frequency, or those arising 
from inflation, that are flat [53]. 



We have plotted in Fig. [TH] the sensitivity of planned 
GW interferometers like LIGO, LISA and BBO, together 
with the present bounds from CMB anisotropics (GUT 
inflation), from Big Bang Nucleosynthesis (BBN) and 
from milisecond pulsars (ms pulsar). Also shown are 
the expected stochastic backgrounds of chaotic inflation 
models like X(p'^, both coupled and pure, as well as the 
predicted background from two different hybrid infla- 
tion models, a high-scale model, with v = 10^'^ Alp and 
X g'^ ^ 0.1, and a low-scale model, with v — 10~''^Afp 
and X ^ ^ 10^^^, corresponding to a rate of expansion 
H ^ 100 GeV. The high-scale hybrid model produces 
typically as much gravitational waves from preheating 
as the chaotic inflation models. The advantage of low- 
scale hybrid models of inflation is that the background 
produced is within reach of future GW detectors like 
BBO [6]. It is speculated that future high frequency laser 
interferometers could be sensitive to a GWB in the MHz 
region jl^], although they are still far from the bound 
marked with an interrogation sign. 
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For a high-scale model of inflation, we may never see 
the predicted GW background coming from preheating, 
in spite of its large amplitude, because it appears at very 
high frequencies, where no detector has yet shown to 
be sufficiently sensitive. On the other hand, if inflation 
occured at low scales, even though we will never have a 
chance to detect the GW produced during inflation in 
the polarization anisotropics of the CMB, we do expect 
gravitational waves from preheating to contribute with 
an important background in sensitive detectors like 
BBO. The detection and characterization of such a GW 
background, coming from the complicated and mostly 
unknown epoch of rehating of the universe, may open a 
new window into the very early universe, while providing 
a new test on inflationary cosmology. 
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